program chem_setstate_test
   use mod_chem
   use iso_fortran_env, only: dp => real64
   implicit none

   type(type_chem_phases) :: phases
   integer, parameter :: n = 1981
   real(dp), allocatable :: T(:), X(:, :), rho_compared(:), H, U
   real(dp) :: p = 101325.0_dp, rho

   integer :: l, i, ll

   integer :: file_io(100)

   phases = chem_file_load("air_5species_NASA9.yaml", "air")
   allocate (T(n), X(5, n), rho_compared(n))
   open (newunit=file_io(100), file="5S_P_xi.plt")
   read (file_io(100), *)
   read (file_io(100), *)
   do l = 1, n
      read (file_io(100), *) T(l), X(2, l), X(5, l), X(1, l), X(3, l), X(4, l), rho_compared(l)
   end do
   close (file_io(100))

   l = 10
   call phases%set_state(T=T(l), X=X(1:5, l), P=p)
   U = phases%get_internal_energy_mass()
   H = phases%get_enthalpy_mass()

   !! test set H
   call phases%set_state(H=H, X=X(1:5, l), P=p)
   print *, T(l), phases%get_temperature()

   !! test set U
   call phases%set_state(U=U, X=X(1:5, l), P=p)
   print *, T(l), phases%get_temperature()

   print *, phases%get_molar_density()

end
